
#include "stdafx.h"
#define MYLIB_API_EXPORT

#include "MeshSimplification.h"
#include "KernelConverter.h"

//initialize the static member
std::auto_ptr<Simplify_func> Simplify_func::Simp_instance(0);

int Simplify_func::Simplification_Turk(ModelOBJ* model,float ratio)
{
	ModelOBJ::const_iterator itr;
	long simNum = 0;
	//for each solid 
	for (itr = model->begin(); itr != model->end(); ++itr)
	{
		Polyhedron_3* cur = (*itr)->Solid_mesh;
		//do not simplify simple models(city50 300)
		float curRatio = ratio;
		// The surface and stop conditions are mandatory arguments.
		// The index maps are needed because the vertices and edges
		// of this surface lack an "id()" field CGAL.
		CGAL::Surface_mesh_simplification::LindstromTurk_params placeparam(0.5,0.5,0);
		CGAL::Surface_mesh_simplification::Count_ratio_stop_predicate<Polyhedron_3> stop(curRatio);
		int r = CGAL::Surface_mesh_simplification::edge_collapse(*cur
			,stop
			,CGAL::vertex_index_map(boost::get(CGAL::vertex_external_index,*cur))
			.edge_index_map (boost::get(CGAL::edge_external_index,*cur))
			.get_cost(CGAL::Surface_mesh_simplification::LindstromTurk_cost<Polyhedron_3>(placeparam))
			.get_placement(CGAL::Surface_mesh_simplification::LindstromTurk_placement<Polyhedron_3>(placeparam))
			);
		//
		if (!cur->is_closed())
		{
			CGAL_assertion(0);
		}
		simNum += r;
	}
	return simNum;
}

int Simplify_func::Simplification_Halfedge(ModelOBJ* model,float ratio)
{
	ModelOBJ::const_iterator itr;
	long simNum = 0;
	//for each solid 
	for (itr = model->begin(); itr != model->end(); ++itr)
	{
		Polyhedron_3* cur = (*itr)->Solid_mesh;
		float curRatio = ratio;
		// The surface and stop conditions are mandatory arguments.
		// The index maps are needed because the vertices and edges
		// of this surface lack an "id()" field CGAL.
		CGAL::Surface_mesh_simplification::LindstromTurk_params placeparam(0.5,0.5,0);
		CGAL::Surface_mesh_simplification::Count_ratio_stop_predicate<Polyhedron_3> stop(curRatio);
		int r = CGAL::Surface_mesh_simplification::edge_collapse(*cur
			,stop
			,CGAL::vertex_index_map(boost::get(CGAL::vertex_external_index,*cur))
			.edge_index_map (boost::get(CGAL::edge_external_index,*cur))
			.get_cost(CGAL::Surface_mesh_simplification::LindstromTurk_cost<Polyhedron_3>(placeparam))
			.get_placement(My_placement<Polyhedron_3>())
			);
		//
		if (!cur->is_closed())
		{
			CGAL_assertion(0);
		}
		simNum += r;
	}
	return simNum;
}


long Simplify_func::eliminate_degenerations(Polyhedron_3& p)
{
	long Num_deg = 0;
	Polyhedron_3::Face_iterator fit;
	for (fit = p.facets_begin(); fit != p.facets_end(); ++fit)
	{
		if (fit->halfedge()->is_border()) continue; 
		//detect generation by facet normal
		Polyhedron_3::Plane_3 plane(fit->halfedge()->vertex()->point(), 
			fit->halfedge()->next()->vertex()->point(), 
			fit->halfedge()->opposite()->vertex()->point() 
			); 
		double a = plane.a();
		double b = plane.b();
		double c = plane.c();
		//handle degeneration(tackle one edge each time)
		//if (plane.is_degenerate())//too strick
		if (fabs(a) < DEGENERATION_TEST && fabs(b) < DEGENERATION_TEST && fabs(c) < DEGENERATION_TEST )
		{
			// equality vertex
			if (fit->halfedge()->vertex()->point() == fit->halfedge()->next()->vertex()->point())
			{
				deg_edge_collapse(p, fit->halfedge()->next());
				fit = p.facets_begin();
				Num_deg++;
				continue;
			}
			else if (fit->halfedge()->opposite()->vertex()->point() == fit->halfedge()->next()->vertex()->point())
			{
				deg_edge_collapse(p, fit->halfedge()->next()->next());
				fit = p.facets_begin();
				Num_deg++;
				continue;
			}
			else if (fit->halfedge()->opposite()->vertex()->point() == fit->halfedge()->vertex()->point())
			{
				deg_edge_collapse(p, fit->halfedge());
				fit = p.facets_begin();
				Num_deg++;
				continue;
			}
			//collinear vertices
			else if (CGAL::collinear(fit->halfedge()->vertex()->point(), 
				fit->halfedge()->next()->vertex()->point(), 
				fit->halfedge()->opposite()->vertex()->point()))
			{
				//delete the shortest edge
				if (CGAL::squared_distance(fit->halfedge()->vertex()->point(),fit->halfedge()->next()->vertex()->point())
					<CGAL::squared_distance(fit->halfedge()->next()->vertex()->point(),fit->halfedge()->opposite()->vertex()->point())
					&& CGAL::squared_distance(fit->halfedge()->vertex()->point(),fit->halfedge()->next()->vertex()->point())
					<CGAL::squared_distance(fit->halfedge()->opposite()->vertex()->point(),fit->halfedge()->vertex()->point()))
				{
					deg_edge_collapse(p, fit->halfedge()->next());
				}
				else if (CGAL::squared_distance(fit->halfedge()->next()->vertex()->point(),fit->halfedge()->opposite()->vertex()->point())
					<CGAL::squared_distance(fit->halfedge()->vertex()->point(),fit->halfedge()->next()->vertex()->point())
					&& CGAL::squared_distance(fit->halfedge()->next()->vertex()->point(),fit->halfedge()->opposite()->vertex()->point())
					<CGAL::squared_distance(fit->halfedge()->opposite()->vertex()->point(),fit->halfedge()->vertex()->point()))
				{
					deg_edge_collapse(p, fit->halfedge()->next()->next());
				} 
				else
				{
					deg_edge_collapse(p, fit->halfedge());
				}

				fit = p.facets_begin();
				Num_deg++;
				continue;
			}
			else
			{
				//try to use tolerance
				// 				double test11 = fit->halfedge()->vertex()->point().x().to_double();
				// 				double test12 = fit->halfedge()->vertex()->point().y().to_double();
				// 				double test13 = fit->halfedge()->vertex()->point().z().to_double();
				// 
				// 				double test21 = fit->halfedge()->opposite()->vertex()->point().x().to_double();
				// 				double test22 = fit->halfedge()->opposite()->vertex()->point().y().to_double();
				// 				double test23 = fit->halfedge()->opposite()->vertex()->point().z().to_double();
				// 
				// 				double test31 = fit->halfedge()->next()->vertex()->point().x().to_double();
				// 				double test32 = fit->halfedge()->next()->vertex()->point().y().to_double();
				// 				double test33 = fit->halfedge()->next()->vertex()->point().z().to_double();
				// 				//
				// 				if (fabs(test11 - test21) < COPLANER_DIF_TOL*1e7 && fabs(test12 - test22) < COPLANER_DIF_TOL*1e7 &&
				// 					fabs(test13 - test23) < COPLANER_DIF_TOL*1e7)
				// 				{
				// 					deg_edge_collapse(p, fit->halfedge());
				// 					fit = p.facets_begin();
				// 					Num_deg++;
				// 					continue;
				// 				}
				// 				else if (fabs(test21 - test31) < COPLANER_DIF_TOL*1e7 && fabs(test22 - test32) < COPLANER_DIF_TOL*1e7 &&
				// 					fabs(test23 - test33) < COPLANER_DIF_TOL*1e7)
				// 				{
				// 					deg_edge_collapse(p, fit->halfedge()->next()->next());
				// 					fit = p.facets_begin();
				// 					Num_deg++;
				// 					continue;
				// 				}
				// 				else if (fabs(test11 - test31) < COPLANER_DIF_TOL*1e7 && fabs(test12 - test32) < COPLANER_DIF_TOL*1e7 &&
				// 					fabs(test13 - test33) < COPLANER_DIF_TOL*1e7)
				// 				{
				// 					deg_edge_collapse(p, fit->halfedge()->next());
				// 					fit = p.facets_begin();
				// 					Num_deg++;
				// 					continue;
				// 				}
				// 				else
				// 				{
				// 					std::cout<<"degeneration could not be handled"<<std::endl;
				// 				}
			}
		}
	} 
	return Num_deg;
}

//should be improved (colinear vertices)
int Simplify_func::merge_coplaner_facets(ModelOBJ* model,double dtol)
{
	ModelOBJ::const_iterator itr;
	long mergeNum = 0;
	long solidNum = 0;
	//for each solid 
	for (itr = model->begin(); itr != model->end(); ++itr)
	{
		Polyhedron_3* cur = (*itr)->Solid_mesh;
		std::cout<<"Solid: "<<solidNum++;
		//
		//eliminate degeneration
		std::cout<<"(Degeneration: "<<eliminate_degenerations(*cur)<<") ";
		//for each vertex
		long preNum = mergeNum;
		Vertex_iterator vit = cur->vertices_begin();
		while(vit != cur->vertices_end())
		{
			if (merge_facets_around_vertex(cur, vit, dtol))
			{
				vit = cur->vertices_begin();
				mergeNum++;
			}
			else
			{
				++vit;
			}
		}
		//collapse vertices
		std::cout<<"(Co-planer Vertice: "<<mergeNum-preNum<<") "<<std::endl;
	}
	return mergeNum;
}

int Simplify_func::clean_nef_polyhedron(Nef_polyhedron_3& poly)
{
	long mergeNum = 0;
	//
	typedef CGAL::Polyhedron_3<EK, My_items> EK_Poly_3;
	EK_Poly_3 ek_polycpy;
	if (poly.is_simple())
	{
		poly.convert_to_polyhedron(ek_polycpy);
	}
	else
	{
		return mergeNum;
	}
	//kernel conversion EK to IK
	
	Polyhedron_3 polycpy;
	if (!CKernelConverter::GetSingleton()->EK_to_IK(&ek_polycpy, &polycpy))
	{
		return 0;
	}
	//
	//remove coplaner vertices
	Vertex_iterator vit = polycpy.vertices_begin();
	while(vit != polycpy.vertices_end())
	{
		if (merge_facets_around_vertex(&polycpy, vit, NORMAL_DIF_TOL))
		{
			vit = polycpy.vertices_begin();
			mergeNum++;
		}
		else
		{
			++vit;
		}
	}
	//eliminate degeneration
	std::cout<<"(Degeneration: "<<eliminate_degenerations(polycpy)<<")";

	//convert back
	if (mergeNum > 0 && polycpy.is_closed())
	{
		//kernel conversion IK to EK
		if (!CKernelConverter::GetSingleton()->IK_to_EK(&polycpy, &ek_polycpy))
		{
			return 0;
		}
		//
		Nef_polyhedron_3 newpoly(ek_polycpy);
		poly = newpoly.regularization();
		return mergeNum;
	}
	else
	{
		return 0;
	}
}

int Simplify_func::clean_polyhedron(Polyhedron_3& poly)
{
	long mergeNum = 0;
	//conversion
	if (!poly.is_closed())
	{
		return mergeNum;
	}
	//remove coplaner vertices
	Vertex_iterator vit = poly.vertices_begin();
	while(vit != poly.vertices_end())
	{
		if (merge_facets_around_vertex(&poly, vit, NORMAL_DIF_TOL))
		{
			vit = poly.vertices_begin();
			mergeNum++;
		}
		else
		{
			++vit;
		}
	}
	//eliminate degeneration
	std::cout<<"(Degeneration: "<<eliminate_degenerations(poly)<<")";
	//
	if (mergeNum > 0 && poly.is_closed())
	{
		return mergeNum;
	}
	else
	{
		return 0;
	}
}
bool Simplify_func::deg_edge_collapse(Polyhedron_3& poly, Halfedge_handle h)
{
	// do the collapse
	if (CGAL::circulator_size(h->vertex_begin()) == 3)
	{
		poly.erase_center_vertex(h);
	}
	else if(CGAL::circulator_size(h->vertex_begin()) > 3)
	{
		//join facets
		Halfedge_handle vth;
		if (CGAL::circulator_size(h->next()->vertex_begin())>=3&&
			CGAL::circulator_size(h->next()->opposite()->vertex_begin())>=3)
		{
			vth = poly.join_facet(h->next());
		}
		else
		{
			std::cout<<"circulator size"/*<<CGAL::circulator_size(vit->halfedge()->next()->vertex_begin())*/<<std::endl;
			return false;
		}
		//
		if (CGAL::circulator_size(vth->opposite()->prev()->vertex_begin())>=3&&
			CGAL::circulator_size(vth->opposite()->prev()->opposite()->vertex_begin())>=3)
		{
			vth = poly.join_facet(vth->opposite()->prev());
		}
		else
		{
			std::cout<<"circulator size"/*<<CGAL::circulator_size(vth->opposite()->prev()->vertex_begin())*/<<std::endl;
			return false;
		}
		//join vertex
		if (CGAL::circulator_size( vth->prev()->facet_begin())>= 4&&
			CGAL::circulator_size( vth->prev()->opposite()->facet_begin())>= 4)
		{
			poly.join_vertex(vth->prev());	
		} 
		else
		{
			std::cout<<"circulator size"/*<<CGAL::circulator_size( vth->prev()->facet_begin())*/<<std::endl;
			return false;
		}

	}
	else
	{
		std::cout<<"halfedge around vertex: "<<CGAL::circulator_size(h->vertex_begin())<<std::endl;
		return false;
	}

	return true;
}

bool Simplify_func::edge_collapse(Polyhedron_3* poly, Halfedge_handle h, const Vector_3& refNormal)
{
	//predict the normals produced by the collapse
	typedef Polyhedron_3::Halfedge_around_vertex_circulator HV_circulator;//clockwise
	HV_circulator hcir = h->vertex_begin();
	int iNum = CGAL::circulator_size(hcir);
	hcir++;//start from the second facet
	for (int i = 0; i < iNum-2; ++i)//predict the other facets
	{
		Polyhedron_3::Plane_3 plane = Polyhedron_3::Plane_3(
			hcir->opposite()->vertex()->point(),
			h->opposite()->vertex()->point(),
			hcir->next()->vertex()->point());
		//
		double dTheta = 1-(plane.orthogonal_vector()*refNormal)/
			sqrt(plane.orthogonal_vector().squared_length()*refNormal.squared_length());
		if ( dTheta > 0.2)
		{
			//collapse will cause flip of normal
			return false;
		}
		// 		if (plane.orthogonal_direction() != refNormal.direction())
		// 		{
		// 			return false;
		// 		}
		++hcir;
	}

	// do the collapse
	if (CGAL::circulator_size(h->vertex_begin()) == 3)
	{
		poly->erase_center_vertex(h);
	}
	else if(CGAL::circulator_size(h->vertex_begin()) > 3)
	{
		//join facets
		Halfedge_handle vth;
		if (CGAL::circulator_size(h->next()->vertex_begin())>=3&&
			CGAL::circulator_size(h->next()->opposite()->vertex_begin())>=3)
		{
			vth = poly->join_facet(h->next());
		}
		else
		{
			std::cout<<"circulator size 1"/*<<CGAL::circulator_size(vit->halfedge()->next()->vertex_begin())*/<<std::endl;
			return false;
		}
		//
		if (CGAL::circulator_size(vth->opposite()->prev()->vertex_begin())>=3&&
			CGAL::circulator_size(vth->opposite()->prev()->opposite()->vertex_begin())>=3)
		{
			vth = poly->join_facet(vth->opposite()->prev());
		}
		else
		{
			std::cout<<"circulator size 2"/*<<CGAL::circulator_size(vth->opposite()->prev()->vertex_begin())*/<<std::endl;
			//undo
			//poly->split_facet(vth, vth->prev()->prev());
			//
			return false;
		}
		//join vertex
		if (CGAL::circulator_size( vth->prev()->facet_begin())>= 4&&
			CGAL::circulator_size( vth->prev()->opposite()->facet_begin())>= 4)
		{
			poly->join_vertex(vth->prev());	
		} 
		else
		{
			std::cout<<"circulator size 3"/*<<CGAL::circulator_size( vth->prev()->facet_begin())*/<<std::endl;
			// 				//undo
			// 				vth = poly->split_facet(vth, vth->next()->next());
			// 				poly->split_facet(vth->next()->opposite(), vth->opposite()->prev()->opposite()->next());
			//
			return false;
		}
	}
	else
	{
		std::cout<<"halfedge around vertex: "<<CGAL::circulator_size(h->vertex_begin())<<std::endl;
		return false;
	}

	return true;
}

bool Simplify_func::merge_facets_around_vertex(Polyhedron_3* poly, Vertex_iterator vit, double dTol)
{
	typedef Polyhedron_3::Halfedge_around_vertex_circulator HV_circulator;
	//
	std::vector<Vector_3> Normal_vec;
	Vector_3 averageNormal(0.0, 0.0, 0.0);
	//for each halfedge
	HV_circulator hcir = vit->vertex_begin();
	do
	{
		//whether the vertex is on the boundary
		if (hcir->is_border()||hcir->opposite()->is_border())
		{
			std::cout<<"Border!!"<<std::endl;
			return false;
		}
		//calculate the normals around the vertex
		Polyhedron_3::Plane_3 plane = Polyhedron_3::Plane_3(hcir->vertex()->point(), 
			hcir->next()->vertex()->point(), 
			hcir->opposite()->vertex()->point());
		Normal_vec.push_back(plane.orthogonal_vector());
		averageNormal = averageNormal + plane.orthogonal_vector();
		++hcir;
	}while(hcir != vit->vertex_begin());
	//compare the maximum normal deviation with the Tol
	averageNormal = averageNormal / (int)Normal_vec.size();
	std::vector<Vector_3>::iterator ditr = Normal_vec.begin();
	double maxDev = (*ditr)*averageNormal/sqrt((*ditr).squared_length()*averageNormal.squared_length());
	while(++ditr != Normal_vec.end())
	{
		double dDev = (*ditr)*averageNormal/sqrt((*ditr).squared_length()*averageNormal.squared_length());
		maxDev = dDev > maxDev ? dDev : maxDev;
	}
	//merge by using edge collapse
	//in order to prevent flip of normal, test all the halfedges around the vertex
	if ((1 - maxDev) < dTol)
	{
		hcir = vit->vertex_begin();//clockwise
		do
		{
			if (edge_collapse(poly, hcir, Normal_vec[0]))
			{
				return true;
			}
			++hcir;
		}while (hcir != vit->vertex_begin());
	} 
	return false;
}